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INTRODUCTION 


The  notion  of  self-organized  critical  systems  has  been  introduced  by  Bak  et  al.  (ref  1)  to 
provide  a  consistent  explanation  for  the  fractal  spatial  structures,  power  law  distributions,  and 
flicker  noise  commonly  observed  in  spatially  extended-driven  dynamical  systems.  They  have 
proposed  a  deposition  model  with  a  geometric  configuration  that  evolves  toward  a  critical  form. 
The  critical  state  is  an  attractor  for  the  dynamics  of  the  system,  and  the  model  demonstrates  how 
complex  fractal  surfaces  evolve  from  simple  system  dynamics. 

APPROACH 

The  cellular  automata  employed  in  this  study  are  based  on  a  two-dimensional  regular 
lattice  of  cells  using  the  rules  given  by  Kadanoff  et  al.  (ref  2).  The  basic  variable  is  z(i,j),  where 
(i,j)  are  the  spatial  indices  and  z  represents  the  height  of  the  lattice.  The  system  is  initialized  to  a 
planar  surface,  z(i,j)  =  0.  Particles  are  added  at  random  locations  z(i,j)  such  that  z(i,j)  —>  z(i,j)  + 

1  as  long  as  the  sites  are  stable.  After  a  particle  is  added  and  after  particle  rearrangements,  the 
stability  of  all  lattice  sites  is  determined  by  z(i,j)  and  the  nearest  neighbor  z  values.  In  the  model 
treated  here,  a  site  becomes  unstable  when 


%  { z(i,  j )  -  z(i,  j  -  k))H  (. zii ,  j )  -  z(i,  j-k))  + 

*=-l 

(z(i,  j )  -  z(i-k,  j))H(z(i ,  j) -  z(i-k, j))}  >  (7C  (1) 

where  H  is  the  Heaviside  function  and  ac  is  the  stability  parameter.  If  site  (i,j)  becomes  unstable, 
z(i,j)  is  reduced  by  Az(i,j),  where  Az(i,j)  = 


%{H(z(i,j)-z(i,j-k))  +  H(z(i,j)-z(i-k,j))}  (2) 

*=— i 

The  height  of  the  corresponding  nearest  neighbor  cells  is  incremented  by  1  if  their  height 
is  less  than  z(i,j).  If  all  sites  are  stable,  a  new  particle  is  added  to  the  system.  The  boundaries  are 
periodic  and  the  local  rule  is  applied  recursively  to  the  cells  whose  state  is  affected  by  the 
unstable  site.  The  diffusion  process  continues  until  there  are  no  more  unstable  sites,  but  no  new 
particles  are  added  until  the  lattice  stabilizes.  The  number  of  sites  that  are  changed  as  the  lattice 
reorganizes  after  a  particle  is  added  is  the  size  of  an  avalanche,  and  the  total  number  of  particles 
that  have  been  added  to  the  system  is  denoted  n. 

The  avalanches  have  been  extensively  explored  (refs  1-4),  and  it  has  been  demonstrated 
that  avalanche  sizes  have  the  form  of  cutoff  hyperbolic  distributions.  However,  quantitative 
analysis  of  the  geometric  scaling  properties  of  the  evolving  automata  structures  has  been  limited. 
One  study  by  Meisel  and  Johnson  (ref  5)  has  explored  the  structures  in  terms  of  self-similar 
fractals,  and  it  has  been  shown  that  there  is  a  monotonic  increase  in  the  Hausdorf-Besokovitch 
dimension  D  with  increasing  ac. 
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However,  the  nature  of  the  deposition  process  and  the  cellular  automata  models  suggest 
that  the  resulting  surface  morphology  belongs  to  a  class  of  fractals  called  surface  fractals  (ref  6), 
which  are  usually  measured  in  terms  of  self-affine  fractal  parameters  (ref  7).  Self-affine  fractals 
are  statistically  invariant  under  anisotropic  dilations,  whereas  self-similar  fractals  are  statistically 
invariant  under  isotropic  dilations  (ref  8).  The  simulations  being  investigated  in  this  study 
pertain  to  surfaces  in  E3  and  are  characterized  by  statistical  invariance  under  transformations  of 
the  form 


{x,y,z}  ->  (Ax, Xy,XHz}  (3) 

where  the  Hurst  exponent,  H,  describes  the  anisotropic  scaling.  In  a  more  general  case,  the  y 
coordinate  would  have  a  coefficient  XK,  and  in  the  case  of  self-similar  fractals,  H  would  have  a 
value  of  1.  The  invariance  expressed  in  equation  (3)  implies  that  any  point  on  a  self-affine 
surface  can  be  represented  in  the  form  {r,  h(r)},  where  the  height  function  h(r)  is  a  single-valued 
function  of  r  =  { x,y}  e  X.  Equation  (3)  applies  over  a  scaling  range  that  is  measured  in  terms  of 
a  parallel  correlation  length,  £//.  The  parallel  correlation  length  is  the  distance  beyond  which 
there  is  no  correlation  in  heights  between  points  on  the  surface.  The  roughness  of  the  surface  (a) 
is  defined  in  terms  of  the  perpendicular  correlation  length,  B,± ,  which  is  related  to  the  rms 

variations  in  h(r).  Assuming  a  homogenous,  self-affine,  surface  structure,  the  height  correlation 
function  is  given  by  (ref  7) 


Ch(r)  =  (  [h(r0 + r )  -  h(r0)f  ^  (4) 

The  Hurst  exponent  H  is  defined  by  the  small  r  variations  in  Ch(r),  i.e., 

H=V2  d(ln( Ch(r))/d(ln(r))(r  «  £/A)  (5) 

where  the  range  of  r  corresponding  to  the  linear  range  in  Ch(r)  is  referred  to  as  the  scaling  range. 
Beyond  the  scaling  range  at  large  r,  the  elevations  become  uncorrelated  and  Ch(r)  is  given  by 

Ch(r)->  2c? (r  »£//)  (6) 

where  cr  is  the  surface  roughness  and  2  a2  is  referred  to  as  the  perpendicular  correlation  length. 
Figure  1  shows  the  double  logarithmic  plot  of  Ch(r)  versus  r  for  a  self-affine  fractal  model. 
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ln(Ch(r)) 


Figure  1.  Ch(r)  for  a  self-affine  fractal  model. 

MEASURING  THE  SCALING  PARAMETERS 

Figures  2a  and  b  show  the  evolving  surface  geometry  of  a  cellular  automaton  model  with 
oc  =  9.  Figure  3  shows  the  height  correlation  function  corresponding  to  the  surface  after  107 
iterations.  Although  the  height  correlation  function  has  the  functional  form  of  a  self-affine 
fractal  model,  estimates  of  £>/  and  H  depend  on  the  range  of  data  used  to  fit  the  linear  region. 


Figure  2a.  Evolution  of  cellular  automaton  model  at  106  iterations(<7c  =  9). 
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Figure  2b.  Evolution  of  cellular  automaton  model  at  3*108  iterations  (crc  =  9) 


Figure  3.  Ch(r)  for  ac  =  9  after  107  iterations 


A  new  systematic  procedure  has  been  developed  that  measures  the  scaling  parameters  by 
fitting  linear  and  polynomial  splines  to  the  data  and  estimating  the  values  of  £y  and  H  using  the 
data  that  minimizes  the  residuals  in  the  fit.  The  linear  range  is  determined  by  fitting  a  straight 
line,/;,  to  the  first  k  points  and  fitting  a  polynomial,^,  to  points  k  through  m.  k  is  the  index  of 
the  data  point  at  which  the  linear  and  curved  segments  meet.  It  is  the  variable  knot  (index)  of  a 
linear/curved  spline  and  is  a  nonlinear  parameter  determined  by  exhaustion,  m  is  the  index  of  the 
data  point  at  which  Ch(r)  approaches  a  constant  The  value  of  m  is  not  critical,  but  should  be 
an  estimate  of  £//  for  computational  efficiency.  The  coefficients  of  the  polynomials  are  evaluated 
with  continuity  enforced  at  k,  and  the  residuals  (e,)  are  computed  over  the  entire  range.  The 
value  of  k  for  which  the  residuals  are  a  minimum  is  selected  as  the  index  for  the  last  point  of  the 
linear  range,  k  is  computed  using  a  quadratic,  cubic,  and  quartic  for/2.  C1  and  C2  continuity  are 
enforced  for  the  quartic  fit.  C2  continuity  is  not  enforced  for  the  cubic  fit  to  allow  some 
overshoot  in  the  data.  The  largest  of  the  three  values  of  k  is  used  to  ensure  the  maximum  number 
of  points  in  the  linear  region  is  selected  to  determine  H.  If  k  is  less  than  3,  it  is  assumed  that  no 
linear  region  exists  and  that  the  system  does  not  adhere  to  fractal  theory.  |//  is  then  determined 
from  the  intersection  of  the  linear  fit  and  The  points  used  to  compute  |j.  are  at  large  r  where 
no  correlation  is  observed.  The  procedure  for  the  quartic  fit  is  outlined  below: 

f(x)  =  a  +  b(x-k)  (7) 

f2(x)  =  c  +  d(x-k)  +  e(x-k)2  +  f(x-k)3  +g(x-k)4  (8) 

For  C2  continuity  at  knot  k 

fi  (*)  =  fi (*0;  fi (k)  =  f2 ( k );  /,' (k)  =  f2  (k) 
and  the  residuals  (£,-)  become 

£i=a+bpi-yi  for  i  <  ik 
=a+bpi+fqi  +  gri-yi  fori>ik 

where 

Pi  -  x,-k;  qt  =  (x,.  -  kf ;  rt  =  (xi  -  k)4  (9) 
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Minimizing  the  sum  of  the  squares  of  the  residuals  S 


1=0  LR  1=0 


dS 

da 


=  0=>a*m  +  b1£pi+f'£q.+g'£r.='£y. 

LR  R  R  LR 


dS 


3r  - 0  =>  a2  Pi + bY  pf + fY  pa + §Y  PiK, = Y  Pi  y> 

0D  LR  LR  R  R  LR 


db 

dS_ 

dg 

dS 


=o=>aYri+bY,  pm  +  fYwi + sY  ri  =  Y  r^i 


—  =  0  =>  a^q,  +b'Y,piqi  +  /J>,2  +  *£$,»}  =YWi 


(10) 

(11) 

(12) 

(13) 

(14) 


The  systems  of  equations  (11)  through  (14)  are  solved  for  coefficients  a,b,f,  and  g  for  k 
ranging  from  2  to  m-1.  The  data  point  corresponding  to  the  value  of  k  for  which  residuals  are  a 
minimum  is  used  as  the  last  point  in  the  linear  fit.  <y/  is  then  determined  from  the  intersection  of 
the  linear  fit  and  In  this  investigation,  k  was  also  computed  using  a  cubic  (g  =  0)  and 
quadratic  (f=g  =  0)  for/2.  The  largest  value  of  k  evaluated  using  the  three  different  polynomials 
for/2  is  used  in  computing  H.  The  following  system  of  equations  is  solved  to  determine  the 
coefficients  of  the  cubic  and  quartic  fits  for/2: 

Cubic  fits: 


fi(x)  =  c  +  d(x-k)  +  e(x-k)2  +  f(x-k)3  (15) 

Pi  =  -  k\  Pi  =  (*,  -  k )2;  r,  =  (x:  - kf  (16) 

dS 

—  =  0  =*a*m  +  b'£pi+ey£qi+fYri=Yyi  (17) 

°a  LR  R  r  LR 

dS 

— =o  =>aYpt+bYpt+eYp&+fYpiri  =!>,*  as> 

0D  LR  LR  R  R  LR 
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Quadratic  fits: 


0  =>  aj^q,  +  b£p,q{  +  e^qf  +  f^qfi  = 


—  =  0  =>  a^r,  + Pirt  +  +  f^r2  = 


f2(x )  =  c  +  d(x  -  k)  +  e(x  -  k)2 


p.=xi-t,qi=(xi-ky 


=  0  =*a*m  +  b2Jpi+e£qi  =  ^yt 


=  0=>a^Pi+bZ  pf  +  pa  =  £  piyi 


LR  LR 


—  =  o  =>  q.  + b£  Pa  +  e X  4?  =  X  ?< 

OZ  R  R  R  R 


Figure  3  shows  the  intermediate  results  with  4  points  (k  =  4)  in  the  linear  region.  The 
final  value  of  k  was  15  using  a  cubic  for/2  which  resulted  in  the  scaling  parameters  of  H  -  0.19, 
and  I// =  8.08  for  ^=8.76. 

DYNAMIC  SCALING  HYPOTHESIS 

Surface  structure  in  deposition  processes  is  time  dependent  (ref  7),  and  a  measure  of 
evolving  morphology  can  be  obtained  by  applying  a  generalization  of  equation  (4)  as 


Ch(r,t)  =  (  [h(r0  +r,t0  +t)~  h(r0 ,  t0)J  ^ 


At  small  t,  due  to  the  absence  of  a  characteristic  time  scale  (ref  7),  £//  and  £1  are  proportional  to 
powers  of  time  t  and 

&(»  oc  t1'!  fort  <  T  (27) 

£>±(1)°*^  for  t  <  T  (28) 
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^  is  a  conventional  exponent  for  finite  size  effects,  and  in  general,  because  of  anisotropy 
in  the  surface,  1  /£  *B.  r  is  the  time  at  which  |/7  is  equal  to  the  size  of  the  sample.  In  this 
investigation,  t  corresponds  to  the  number  of  particles  added  to  the  system,  n. 

MPI  IMPLEMENTATION 

The  generalized  height  correlation  function,  Ch(r,t),  provides  useful  information  about  the 
dynamics  of  growth  processes,  but  is  numerically  expensive.  Therefore,  a  suite  of  homogenous 
parallel  procedures  has  been  developed  to  map  the  solution  to  a  multicomputer  platform  using  an 
implementation  of  the  Message  Passing  Interface  (MPI)  called  LAM  (ref  9)  for  interprocess 
communication.  LAM  is  a  full  implementation  of  the  MPI  standard,  supports  heterogeneous 
computer  networks,  and  provides  direct  communication  between  application  processes.  It  is  a 
node-oriented  computing  environment  that  uses  a  unique  identifier  assigned  to  each  node  as  the 
primary  synchronization  mechanism.  The  nodes  are  usually  fully  connected  (maximum  1-hop 
distance)  since  the  network  is  generally  a  shared  resource.  In  our  implementation,  a  master 
process  assigns  equal-sized  regions  of  the  lattice  to  slave  processing  nodes  for  computation  of 
partial  Crfr,t)  data.  The  processing  nodes  return  the  results  to  the  master  process  where  Ch{r,t)  is 
updated  and  the  slave  processes  assigned  new  regions  to  analyze.  Although  this  approach  does 
not  exploit  the  communications  advantages  of  a  fully  connected  network,  the  load  balancing 
results  in  a  performance  improvement  directly  proportional  to  the  number  of  workstations.  This 
is  likely  due  to  the  high  computation/communication  ratio  and  the  relatively  small  number  of 
workstations  making  up  the  multicomputer. 

RESULTS 

The  dynamics  of  the  evolving  surface  morphology  of  the  cellular  automata  models  with 
<yc  =  9,  10,  11,  13,  and  15  on  a  256x256  lattice  were  characterized  in  terms  of  the  spatial  scaling 
parameters,  H,  rf/,  and  Q± .  The  height  correlation  function  was  computed  at  discrete  intervals 
(iterations),  and  the  largest  value  of  k  resulting  in  the  minimum  residuals  using  three  different 
polynomials  for/2  was  computed.  The  spatial  scaling  measures  were  determined  for  those 
surfaces  for  which  fractal  theory  applied  (k  >2).  These  values  were  used  to  determine  the 
dynamic  scaling  exponents  1/C  and  B.  Table  1  shows  the  measured  values  of#,  £// ,  and  £*  for 
trc  =  9  at  each  interval  and  the  degree  of  the  polynomial  f2  that  resulted  in  the  largest  value  of  it. 
The  table  shows  that  H  remains  relatively  constant  as  the  surface  evolves.  Results  are  similar  for 
all  other  <7C.  In  all  cases  studied,  a  quartic  fit  for  f2  resulted  in  the  largest  linear  region  only  once. 


Table  1.  Spatial  Scaling  Parameters  with  oc  =  9 


k 

f2 

H 

mm 

mm 

i 

3 

3 

0.16 

3.13 

5.42 

3 

7 

2 

0.17 

5.37 

6.75 

5 

12 

3 

0.17 

6.69 

7.46 

7 

13 

3 

0.18 

7.32 

8.08 

10 

15 

3 

0.19 

8.08 

8.76 

30 

19 

3 

0.21 

11.59 

10.80 

50 

0.22 

13.20 

12.06 

70 

26 

2 

0.22 

100 

37 

3 

0.22 

15.33 

12.81 

300 

35 

2 

0.23 

24.78 

16.61 

The  evolving  surface  morphologies  of  the  cellular  automata  models  adhere  well  to  the 
dynamic  scaling  hypothesis  as  shown  in  Figure  4.  The  figure  shows  typical  In-ln  plots  of  spatial 
scaling  parameters  from  which  the  dynamic  scaling  exponents  are  derived.  The  exponents 
l/£  and  B  are  determined  from  the  slope  of  the  straight-line  fit  using  all  of  the  data  points.  The 
results  for  this  and  the  other  models  are  summarized  in  Table  2.  The  data  show  that  the  value  of 
H  increases  with  increasing  crc  and  at  any  given  point  in  time,  the  parallel  correlation  length 
decreases  with  increasing  <7C.  This  result  is  expected  because  as  <xc  increases,  more  particles  are 
required  for  an  avalanche  to  occur.  There  is  no  evident  correlation  between  B  and  crc. 


Figure  4.  Parallel  correlation  of  length  versus  time. 
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Table  2.  Dynamic  Scaling  Parameters 


Oc 

_ <«> 

y; 

B 

9 

0.20  ±0.03 

0.33 

0.19 

10 

0.25  ±0.04 

0.29 

0.23 

11 

r  0.30  ±0.03 

0.27 

0.23 

13 

0.36  ±0.03 

0.25 

0.24 

15 

0.41  ±0.01 

0.21 

0.20 

SUMMARY 

The  deposition  of  adherent,  erosion  and  corrosion-resistant  coatings  is  critical  to  the 
performance  of  many  systems.  In  order  to  develop  optimal  coatings,  it  is  important  to 
understand  the  effects  of  changing  the  parameters  controlling  the  deposition/sputtering 
processes.  Models  of  the  deposition  processes  and  quantitative  measures  of  the  condition  of  the 
substrate  surface  and  evolving  coating  are,  at  present,  extremely  limited.  In  this  study,  cellular 
automata  simulations  were  used  to  investigate  the  nature  of  the  vapor  deposition  process  by 
exploring  the  natural  evolution  of  dynamical  dissipative  systems  through  self-organized  critical 
system  analysis  and  spatial  scaling  measures.  A  new  numerical  technique  was  developed  to 
analyze  the  intrinsic  structure  of  evolving  surface  morphology  in  an  effort  to  better  understand 
the  dynamics  of  the  growth  processes.  The  algorithm  is  numerically  expensive,  therefore  the 
computational  problem  has  been  mapped  to  a  multicomputer  platform  using  an  implementation 
of  MPI  call  LAM  for  interprocess  communication.  Using  this  approach,  the  deposition  models 
were  shown  to  agree  well  with  the  dynamic  scaling  theory.  The  technique  is  also  being  used  to 
validate  the  integrity  of  other  deposition  models  through  a  comparative  analysis  with 
experimental  data,  and  to  determine  if  a  correlation  exists  between  intrinsic  surface  structure, 
evolving  surface  morphology,  and  parameters  controlling  the  deposition  process. 
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